source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/SLE_TRPM2_MfMo/'

gene_of_int <- c('TRPM2','HMGB2','CCR2','CD14','TNF')

rey_raw <-
data.table::fread('mission/SLE_TRPM2_MfMo/data/reyne20sepsis/scp_gex_matrix_raw.csv.gz')

rey_raw |> nose()

rey_meta <-
  read_delim('mission/SLE_TRPM2_MfMo/data/reyne20sepsis/scp_meta_updated.txt')

rey_meta <- rey_meta |>
  filter(Cell_Type != 'group') |>
  column_to_rownames('NAME ')

sobj <- rey_raw |>
  column_to_rownames('GENE') |>
  CreateSeuratObject(meta.data = rey_meta, min.cells = 3, min.features = 200)

sobj <- sobj |>
  PercentageFeatureSet('^MT-', col.name = 'mito_ratio')

sobj |> VlnPlot('mito_ratio', group.by = 'Cohort', pt.size = 0)

sobj <- sobj |>
  filter(mito_ratio < 10)

sobj <- sobj |>
  NormalizeData()

sobj <- sobj |>
  filter(Cell_Type != 'Megakaryocyte', Cohort != 'ICU-NoSEP')

sobj <- sobj |>
  mutate(group = case_match(Cohort,
                            'Control' ~ 'HC',
                            c('Leuk-UTI') ~ 'Infection',
                            .default = 'Sepsis'))

# TRPM2 dotplot ---------
cd14mo_grp <- sobj_mo |>
  filter(manual_fine == 'CD14+ Mono') |>
  colnames()

sobj <- sobj |>
  mutate(manual_fine = case_when(.cell %in% cd14mo_grp ~ 'CD14+ Mono',
                                 Cell_Type == 'Mono' ~ 'CD16+ Mono',
                                 .default = Cell_Type))

sobj |>
  DotPlot2d('TRPM2', Cell_Type, group) +
  labs(x = 'Cell type', y = 'Cohort',
       title = 'PBMC TRPM2')

sobj |>
  DotPlot2d('TRPM2', manual_fine, group) +
  labs(x = 'Cell type', y = 'Cohort',
       title = 'PBMC TRPM2') +
  RotatedAxis()

sobj |>
  DotPlot2d('TRPM2', manual_fine, Cohort) +
  labs(x = 'Cell type', y = 'Cohort',
       title = 'PBMC TRPM2') +
  RotatedAxis()

# big tsne ------------
tsne <-
  read_delim('mission/SLE_TRPM2_MfMo/data/reyne20sepsis/scp_tsne_updated.txt',
             col_types = 'cdd') |> na.omit()

sobj <-
  read_rds('mission/SLE_TRPM2_MfMo/data/reyne20sepsis/pbmc.rds')

sobj <- tsne |>
  mutate(.cell = NAME, tSNE1 = X, tSNE2 = Y, .keep = 'none') |>
  left_join(x = sobj, y = _)

sobj |>
  ggplot(aes(tSNE1, tSNE2, fill = Cell_Type)) +
  geom_bin_2d(bins = 256) +
  theme_pubr() +
  labs(title = 'Sepsis & UTI patient PBMC') +
  facet_wrap(~group)

sobj |>
  write_rds('mission/SLE_TRPM2_MfMo/data/reyne20sepsis/pbmc.rds')

# monocyte subset ------------
sobj_mo <- sobj |>
  filter(Cell_Type == 'Mono')

sobj_mo$orig.ident <- sobj_mo$donor_id

## gene of int logfc in total mono ------------
goi_totalmo <- sobj_mo |>
  DotPlot(gene_of_int, group.by = 'donor_id') |>
  pluck('data') |>
  as_tibble() |>
  mutate(donor_id = as.character(id)) |>
  left_join(cohort_meta) |>
  mutate(group = case_match(Cohort,
                            'Control' ~ 'HC',
                            c('Leuk-UTI') ~ 'Infection',
                            .default = 'Sepsis'),
         Cohort = fct_relevel(Cohort, 'Control', 'Leuk-UTI', 'Int-URO', 'URO'))

goi_totalmo |>
  ggplot(aes(Cohort, avg.exp, fill = Cohort)) +
  stat_mean(geom = 'col') +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~features.plot, scales = 'free_y') +
  theme_pubr(x.text.angle = 90) +
  scale_fill_viridis_d(option = 'turbo', begin = .1, end = .9) +
  labs(y = 'Average expression', title = 'Total monocyte in PBMC') +
  stat_compare_means(ref.group = 'Control', label = 'p.signif') +
  scale_y_continuous(expand = expansion(mult = c(NA, .1)))

goi_submo <- sobj_mo |>
  DotPlot2d(gene_of_int, donor_id, manual_fine) |>
  pluck('data') |>
  as_tibble() |>
  mutate(donor_id = group.x) |>
  left_join(cohort_meta) |>
  mutate(group = case_match(Cohort,
                            'Control' ~ 'HC',
                            c('Leuk-UTI') ~ 'Infection',
                            .default = 'Sepsis'),
         Cohort = fct_relevel(Cohort, 'Control', 'Leuk-UTI', 'Int-URO', 'URO'))

goi_submo |>
  filter(group.y == 'CD14+ Mono') |>
  ggplot(aes(Cohort, avg.exp, fill = Cohort)) +
  stat_mean(geom = 'col') +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~features.plot, scales = 'free_y') +
  theme_pubr(x.text.angle = 45) +
  scale_fill_viridis_d(option = 'turbo', begin = .1, end = .9) +
  labs(y = 'Average expression', title = 'CD14+ monocyte in PBMC') +
  stat_compare_means(ref.group = 'Control', label = 'p.signif', vjust = -.5) +
  scale_y_continuous(expand = expansion(mult = c(NA, .3)))

goi_submo |>
  filter(group.y != 'CD14+ Mono') |>
  ggplot(aes(Cohort, avg.exp, fill = Cohort)) +
  stat_mean(geom = 'col') +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~features.plot, scales = 'free_y') +
  theme_pubr(x.text.angle = 45) +
  scale_fill_viridis_d(option = 'turbo', begin = .1, end = .9) +
  labs(y = 'Average expression', title = 'CD16+ monocyte in PBMC') +
  stat_compare_means(ref.group = 'Control', label = 'p.signif', vjust = -.5) +
  scale_y_continuous(expand = expansion(mult = c(NA, .3)))

sobj_mo <- sobj_mo |>
  quick_process_seurat(skip_norm = T)

sobj_mo <- sobj_mo |>
  filter(seurat_clusters != 13)

m2cd14 <- sobj_mo |>
  DotPlot(gene_of_int, cols = 'RdBu') |>
  pluck('data')

m2cd14 |>
  ggplot(aes(features.plot, id, fill = avg.exp.scaled)) +
  geom_tile(color = 'black') +
  scale_fill_distiller(palette = 'RdBu') +
  theme_minimal() +
  labs(title = 'Gene expression in total monocytes',
       y = 'Monocyte clusters', x = 'Gene', fill = 'Scaled expression')

m2cd14 |>
  filter(id != 13) |>
  pivot_wider(names_from = features.plot, values_from = avg.exp, id_cols = id) |>
  ggplot(aes(CD14, TRPM2)) +
  stat_smooth(method = 'lm', color = 'grey', linetype = 'dashed') +
  geom_point(aes(color = CD14 > 4)) +
  geom_text_repel(aes(label = id)) +
  stat_cor() +
  theme_bw() +
  scale_color_hue(labels = c('CD14- Mono', 'CD14+ Mono'), direction = -1) +
  labs(title = 'TRPM2-CD14 correlation in sepsis PBMC monocytes',
       color = 'Subtype')

sobj_mo <- sobj_mo |>
  mutate(manual_fine = ifelse(seurat_clusters %in% c(10,7,5,12,8),
                              'CD14- Mono', 'CD14+ Mono'))

sobj_mo |>
  write_rds('mission/SLE_TRPM2_MfMo/data/reyne20sepsis/pbmc_mono.rds')

sobj_mo <-
  read_rds('mission/SLE_TRPM2_MfMo/data/reyne20sepsis/pbmc_mono.rds')

## frac change -----------
cohort_meta <- sobj |>
  distinct(Cohort, donor_id) |>
  as_tibble()

cd14m_bysample <- sobj |>
  filter(biosample_id == 'CD45') |>
  calc_frac_conf_on_grouped_count(donor_id, manual_fine)

cd14m_bysample |>
  left_join(cohort_meta) |>
  filter(str_detect(manual_fine, 'Mono')) |>
  mutate(group = case_match(Cohort,
                            'Control' ~ 'HC',
                            c('Leuk-UTI') ~ 'Infection',
                            .default = 'Sepsis'),
         Cohort = fct_relevel(Cohort, 'Control', 'Leuk-UTI', 'Int-URO', 'URO')) |>
  ggplot(aes(Cohort, fraction*100, fill = Cohort)) +
  stat_mean(geom = 'col') +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~manual_fine) +
  theme_pubr(x.text.angle = 45) +
  scale_fill_viridis_d(option = 'turbo', begin = .1, end = .9) +
  labs(y = '% in PBMC', title = 'Monocyte subtype proportion in PBMC') +
  stat_compare_means(ref.group = 'Control', label = 'p.signif') +
  scale_y_continuous(expand = expansion(mult = c(NA, .15)))

# COMBAT22 dataset ------------
comb_obs <-
  data.table::fread('mission/SLE_TRPM2_MfMo/data/combat22sepsis/combat_meta.csv.gz')

sep_hc_bool <- comb_obs[, Source %in% c('HV','Sepsis')]

adata <-
  anndata::read_h5ad('mission/SLE_TRPM2_MfMo/data/combat22sepsis/COMBAT-CITESeq-DATA.h5ad')

combat_mex <- adata$layers[['raw']]

combat_mex |> glimpse()

sep_mex <- combat_mex[sep_hc_bool,]

sep_mex <- sep_mex |> t()

sep_mex |> glimpse()

sep_mex |>
  DropletUtils::write10xCounts(
    path = 'mission/SLE_TRPM2_MfMo/data/combat22sepsis/combat_sepsis.h5',x = _)

sep_mex <- Read10X_h5('mission/SLE_TRPM2_MfMo/data/combat22sepsis/combat_sepsis.h5')

sep_mex |> glimpse()

sepsis_meta <- comb_obs |>
  as_tibble() |>
  select(V1, Annotation_cell_type, Annotation_major_subset, QC_scrub_doublet_scores,
         Annotation_minor_subset, GEX_region, scRNASeq_sample_ID, Source) |>
  filter(Source %in% c('HV','Sepsis')) |>
  column_to_rownames('V1')

comb_obs |>
  ggplot(aes(QC_scrub_doublet_scores)) +
  geom_density()

comb_obs |>
  ggplot(aes(GEX_region, QC_scrub_doublet_scores)) +
  geom_boxplot()

sobj <- sep_mex |>
  CreateSeuratObject(names.field = 3, names.delim = '-', min.cells = 3,
                     min.features = 200, meta.data = sepsis_meta)

sobj <- sobj |>
  PercentageFeatureSet('^MT-', col.name = 'mito_ratio')

sobj |> VlnPlot('mito_ratio', pt.size = 0)

sobj <- sobj |> 
  filter(GEX_region != 'E: Doublets', QC_scrub_doublet_scores < .2) |>
  NormalizeData()

combat_umap <- read_csv('mission/SLE_TRPM2_MfMo/data/combat22sepsis/umap.csv')

combat_umap$.cell <- comb_obs$V1

combat_umap <- combat_umap |> select(.cell, X_umap1, X_umap2)

sobj <- sobj |>
  left_join(combat_umap)

## save full rds --------
sobj |>
  write_zstd_rds('mission/SLE_TRPM2_MfMo/data/combat22sepsis/cmbt22sep_zst.rds')

sobj <- 
  read_zstd_rds('mission/SLE_TRPM2_MfMo/data/combat22sepsis/cmbt22sep_zst.rds')

## big umap ---------
sobj |>
  ggplot(aes(X_umap1, X_umap2, fill = cell_type)) +
  geom_bin2d(bins = 512) +
  scale_fill_brewer(palette = 'Paired')

## TRPM2 dotplot -------
sobj |>
  DotPlot2d('TRPM2', Source, major_subset) +
  labs(x = 'Group', y = 'Cell type', title = 'TRPM2 in PBMC')

sobj <- sobj |>
  filter(cell_type != 'nan', major_subset != 'nan', minor_subset != 'nan')

sobj <- sobj |>
  mutate(manual_main = case_when(
    major_subset %in% c('cMono','ncMono','DC') ~ major_subset,
    cell_type == 'PLT' ~ 'Platlet',
    cell_type == 'PB' ~ 'Plasmablast',
    cell_type == 'ERYTH' ~ 'RBC',
    .default = cell_type))

sobj |>
  DotPlot2d('TRPM2', Source, manual_main) +
  scale_radius(range = c(1,3)) +
  labs(x = 'Group', y = 'Cell type', title = 'TRPM2 in PBMC') +
  theme_jpub()

publish_source_plot('sepsis_combat22_pbmc_m2_dotplot')

sobj <- sobj |>
  filter(!(manual_main %in% c('RBC','Platlet','HSC')))

# 24min for 234k cell
sobj <- sobj |>
  quick_process_seurat(skip_norm = T, leiden = F)

rownames(sobj) |> str_subset('^AB-')

sobj |>
  ggplot(aes(umap_1, umap_2, fill = manual_main)) +
  geom_bin2d(bins = 512) +
  theme_jpub(theme_classic) +
  labs(title = 'Sepsis PBMC', fill = 'Cell type') +
  scale_fill_brewer(palette = 'Paired')

publish_pdf('sepsis_combat22_pbmc_umap.pdf', width = 60)

frac_by_sample <- sobj |>
  calc_frac_conf_on_grouped_count(scRNASeq_sample_ID, manual_main)

frac_by_sample <- sobj_mo |>
  distinct(scRNASeq_sample_ID, Source) |>
  right_join(frac_by_sample)

frac_by_sample |>
  mutate(Group = ifelse(Source == 'HV', 'HC', 'Sepsis')) |>
  filter(str_detect(manual_main, 'Mono')) |>
  ggplot(aes(Group, fraction*100, color = Group)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  stat_compare_means(comparisons = list(c('HC','Sepsis')), method = 't.test',
                     label = 'p.signif') +
  scale_color_hue(direction = -1) +
  facet_wrap(~manual_main, scales = 'free_y') +
  scale_y_continuous(expand = expansion(mult = c(NA,.1))) +
  labs(title = 'Cell type proportion', y = '% in CD45+ PBMC') +
  theme_jpub()

publish_source_plot('sepsis_combat22_mono_sub_frac_in_pbmc', width = 100)

## mono subset -------
sobj_mo <- sobj |>
  filter(str_detect(manual_main, 'Mono'))

sobj_mo |>
  write_zstd_rds('mission/SLE_TRPM2_MfMo/data/combat22sepsis_mo.zst.rds')

sobj_mo <-
  read_zstd_rds('mission/SLE_TRPM2_MfMo/data/combat22sepsis_mo.zst.rds')

goi_bysample <- sobj_mo |>
  DotPlot(gene_of_int, group.by = 'scRNASeq_sample_ID') |>
  pluck('data') |>
  as_tibble()

goi_bysample <- sobj_mo |>
  distinct(scRNASeq_sample_ID, Source) |>
  mutate(id = fct(scRNASeq_sample_ID)) |>
  right_join(goi_bysample)

### total mono gene sepsis vs HC ------------
goi_bysample |>
  mutate(Group = ifelse(Source == 'HV', 'HC', 'Sepsis')) |>
  ggplot(aes(Group, avg.exp, color = Group)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  stat_compare_means(comparisons = list(c('HC','Sepsis')), method = 't.test',
                     label = 'p.signif', vjust = .5) +
  scale_color_hue(direction = -1) +
  facet_wrap(~features.plot, scales = 'free_y') +
  scale_y_continuous(expand = expansion(mult = c(NA,.1))) +
  labs(title = 'Total monocytes expression', y = 'Average expression') +
  theme_jpub()

publish_source_plot('sepsis_combat22_totalmono_gene_boxplot',
                    width = 100, height = 80)

sobj_mo <- sobj_mo |>
  quick_process_seurat(skip_norm = T)

inflam15 <- read_csv('mission/SLE_TRPM2_MfMo/results/top15.inflamm.gene.csv')

sobj_mo <- sobj_mo |>
  AddModuleScore(features = list(inflam15$gene), name = 'inflam')

sobj_mo |>
  DotPlot(c(gene_of_int, 'inflam1'), cols = 'RdBu', cluster.idents = T)

goi_bymoclus <- last_plot() |>
  pluck('data') |>
  as_tibble()

goi_bymoclus |>
  BubblePlot() +
  scale_radius(range = c(1,2.5)) +
  labs(title = 'Monocyte clusters') +
  theme_jpub() +
  RotatedAxis()

publish_source_plot('sepsis_combat22_mono_trpm2_dotplot')

### TRPM2 x CCR2/CD14 corr -----------
goi_bymoclus |>
  pivot_wider(id_cols = id, names_from = features.plot, values_from = avg.exp) |>
  mutate(subset = ifelse(id %in% c(4,2,10,15,9,5,8,16,11), 'CD14+ Mono', 'CD16+ Mono')) |>
  ggplot(aes(TRPM2, CCR2)) +
  geom_smooth(method = 'lm', color = 'grey', linetype = 'dashed') +
  geom_point(aes(color = subset)) +
  geom_text_repel(aes(label = id), size = 2) +
  stat_cor(output.type = 'tex', size = 2, label.y = .75) +
  labs(title = 'TRPM2-CCR2 expression correlation in monocytes') +
  theme_jpub()

publish_source_plot('sepsis_combat22_mono_m2_ccr2_cor', width = 70)

goi_bymoclus |>
  pivot_wider(id_cols = id, names_from = features.plot, values_from = avg.exp) |>
  mutate(subset = ifelse(id %in% c(4,2,10,15,9,5,8,16,11), 'CD14+ Mono', 'CD16+ Mono')) |>
  ggplot(aes(TRPM2, CD14)) +
  geom_smooth(method = 'lm', color = 'grey', linetype = 'dashed') +
  geom_point(aes(color = subset)) +
  geom_text_repel(aes(label = id), size = 2) +
  stat_cor(output.type = 'tex', size = 2, label.y = 8) +
  labs(title = 'TRPM2-CD14 expression correlation in monocytes') +
  theme_jpub()

publish_source_plot('sepsis_combat22_mono_m2_cd14_cor', width = 70)

goi_bymoclus |>
  pivot_wider(id_cols = id, names_from = features.plot, values_from = avg.exp) |>
  mutate(subset = ifelse(id %in% c(4,2,10,15,9,5,8,16,11), 'CD14+ Mono', 'CD16+ Mono')) |>
  ggplot(aes(TRPM2, inflam1)) +
  geom_smooth(method = 'lm', color = 'grey', linetype = 'dashed') +
  geom_point(aes(color = subset)) +
  geom_text_repel(aes(label = id), size = 2) +
  stat_cor(output.type = 'tex', size = 2, label.y = 1.8) +
  labs(title = 'TRPM2-CD14 expression correlation in monocytes',
       y = 'Inflammatory score') +
  theme_jpub()

publish_source_plot('sepsis_combat22_mono_m2_inflam_cor', width = 70)

### mono umap --------
sobj_mo |>
  ggplot(aes(umap_1, umap_2, fill = seurat_clusters)) +
  geom_bin2d(bins = 512) +
  scale_fill_manual(values = DiscretePalette(36)) +
  labs(title = 'Monocyte clusters', fill = 'Clusters') +
  theme_jpub(theme_classic)

sobj_mo |>
  DimPlot(cols = DiscretePalette(36), label = T, label.box = T, label.size = 2) +
  ggtitle('Monocytes') +
  theme_jpub(theme_classic) +
  NoLegend()

publish_pdf('sepsis_combat22_mono_leiden_umap.pdf')

sobj_mo |>
  ggplot(aes(umap_1, umap_2, fill = manual_main)) +
  geom_bin2d(bins = 512) +
  labs(title = 'Monocyte clusters', fill = 'Subset') +
  theme_jpub(theme_classic)

publish_pdf('sepsis_combat22_mono_subset_umap.pdf')

sobj_mo |>
  get_abundance_sc_wide('TRPM2') |>
  left_join(x = sobj_mo, y = _) |>
  ggplot(aes(umap_1, umap_2, weight = TRPM2)) +
  geom_bin2d(bins = 256) +
  scale_fill_gradient(low = 'lightgrey', high = 'red') +
  facet_wrap(~Source) +
  theme_jpub(theme_classic) +
  labs(title = 'TRPM2', fill = 'Expression')

publish_pdf('sepsis_combat22_mono_trpm2_featureplot.pdf', width = 100)

### GSEA ----------
cmono_vs_nc_deg <- sobj_mo |>
  FindMarkersAcrossVar(split.by = 'Source', group.by = 'manual_main',
                       ident.1 = 'cMono')

cmono_vs_nc_deg |>
  write_source_csv('sepsis_combat22_cmono_vs_nc_deg')

cmono_vs_nc_deg <- read_csv('mission/SLE_TRPM2_MfMo/results/sepsis_combat22_cmono_vs_nc_deg.csv')

sepsis_cvsnc_gsego <- cmono_vs_nc_deg |>
  filter(cluster == 'Sepsis', p_val_adj < .05) |>
  batch_enrich_path(path = 'GO', method = 'GSEA')

sepsis_cvsnc_gsegosim <- sepsis_cvsnc_gsego[[1]] |>
  simplify()

set.seed(42)

sepsis_cvsnc_gsego[[1]]@result |>
  filter(NES > 0, ONTOLOGY == 'BP') |>
  plot_enrichment(force_regex = 'chronic inflamm')

last_plot() +
  theme_jpub() +
  labs(title = 'GO BP GSEA enriched pathways in sepsis\nCD14+ monocytes vs CD16+ monocytes',
       y = 'Pathway')

publish_source_plot('sepsis_combat22_sepsis_cMono_vs_ncMono_gsea', width = 70)

sepsis_cvsnc_gsego <-
  read_csv('mission/SLE_TRPM2_MfMo/results/sepsis_combat22_sepsis_cMono_vs_ncMono_gsea.csv')

### inflammatory violin ---------
sepsis_infl_gene <- sepsis_cvsnc_gsego |>
  filter(str_detect(Description, 'chronic inflam')) |>
  pull(core_enrichment) |>
  str_split_1('/')

cmono_vs_nc_deg <-
  read_csv('mission/SLE_TRPM2_MfMo/results/sepsis_combat22_cmono_vs_nc_deg.csv')

inflam_go_gene <- map_go_gene('GO:0006954')

sepsis_infl15 <- cmono_vs_nc_deg |>
  filter(cluster == 'Sepsis', avg_log2FC > 0, p_val_adj < .05, gene %in% inflam_go_gene$SYMBOL) |>
  slice_min(p_val_adj, n = 15)

select <- dplyr::select

sobj_mo |>
  filter(Source == 'Sepsis') |>
  bill.violin(sepsis_infl15$gene, manual_main, facet.ncol = 5, pubsize = T) +
  labs(x = 'Cell type', y = 'Normalized expression', fill = 'Cell type',
       title = 'Inflammatory response genes in sepsis monocytes')

publish_pdf('sepsis_combat22_sepsis_cMono_vs_ncMono_inflamm_violin.pdf', width = 90)
